CD Schooolhouse Version 10.0/CD Schooolhouse Version 10.0.ISO/pc/dos/misc/esptest/espprob.pas

Pascal/Delphi Source File | 1994-05-20 | 12.2 KB | 347 lines

UNIT ESPPROB;
{**********************************************************}
{**                                                      **}
{**      A Unit to compute probability  for ESPTEST      **}
{**      Copyright 1991 Phil Mosier                      **}
{**                                                      **}
{**      Used with permision                             **}
{**      From BIO-STATISTICAL MICROCOMPUTING IN PASCAL   **}
{**      BY Klotz and Meyer                              **}
{**      (c) 1985 Rowman & Allanheld Publishers          **}
{**                                                      **}
{**********************************************************}
Interface
Uses Graph, Crt;
Procedure BINOMIAL(Var PROB_B: Real;
                       TRIALS: Word;
                       HITS:   Word);

Procedure CHISQUARE(Var PROB_C: Real;
                        CHI_SQ: Real;
                        DEGREE_FREEDOM: Real);

Implementation
Procedure BINOMIAL(Var PROB_B: Real;
                       TRIALS: Word;
                       HITS:   Word);

{**********************************************************}
{**    From BIO-STATISTICAL MICROCOMPUTING IN PASCAL     **}
{**    BY Klotz and Meyer                                **}
{**********************************************************}

Var
X,XL,A,B,N,P,PL : Real;
{*********************************************************}
{**                                                     **}
{**                                                     **}
{*********************************************************}
Function BINET(A:Real): Real; { Benet's Function }
   Var
   S: Real;
   Begin {BINET}
   {*************************  BINET  ********************}
      If A < 1.5 Then
         BINET := 8.1061467E-2
        {*************************************************}
        {**      Case Binet(1) = 1.5*ln(2pi)            **}
        {*************************************************}
     Else
        {*************************************************}
        {**    Continued Fraction                       **}
        {*************************************************}
        Begin
           S := A + 1.5174736 /( A + 2.269489);
           S := A + ( 195 /371 ) /( A + 22999 /22737 /S);
           BINET := 1 /12 /( A + ( 1 /30 ) /( A + ( 52 /210 ) /S));
        End;
End; {BINET}

Function G(U : Real) : Real;
{*********************************************************}
{**   LN(1+U)-U With Accuracy                           **}
{*********************************************************}
Var
R: Real;
Begin {G}
   If ((U < -0.5) or ( U > 1.0)) Then
      G := Ln(1 + U) - U
   Else
   {*************************************************}
   {**   continued fraction                        **}
   {*************************************************}
   Begin
      R:= 6 + 16 * U /(7 + 9 * U /(8 + 25 * U /(9 + 16 * U /10)));
      G := -U * U /(2 + 4 * U /(3 + U /(4 + 9 * U /(5 +4 * U /R ))));
   End;
End;  { G}

Function ICBETA(X,XL,A,B : Real) : Real;
{*********************************************************}
{**  Program for the incomplete beta distribution.      **}
{**  The algorithm of H.E. Soper (1921), Tracts for     **}
{**  Computer No. VII, Karl Person Ed., Cambridge Univ. **}
{** Press, is used with a,b incremented if either < 1.  **}
{*********************************************************}
Var
C,FXAB : Real;
   Function F(X,XL,A,B : Real) : Real;
   Var
   S,U,V,GU,GV : Real;
      Begin {F}
      U := X * (A + B) /A - 1;
      V := XL * A /B - X;
      If U < -0.9 Then GU := Ln(X * (A + B) /A) - U
      Else GU := G(U);
      If V < -0.9 Then GV := Ln(XL * (A + B) /B) - V
      Else GV := G(V);
      S := 0.5 * Ln(A * B /(( A + B) * 6.2831853));
      S := S + A * GU + B * GV + BINET(A + B) - BINET(A) - BINET(B) - Ln(A);
      If S < - 85.19
         {************************************************}
         {**  Ln(IE-37)                                 **}
         {************************************************}
         Then F := 0
         Else F := Exp(S);
   End; {F}

   Function SOPER(X,XL,A,B : Real) : Real;
   {********************************************************}
   {**                                                    **}
   {********************************************************}
   Const
   EPS = 1.0E-7;
   Var
   V,AI,CK,SUM,ABL,AS,A12,B2,AB4,R : Real;
   S,SMAX,I,K : Integer;
   Begin { SOPER  }
   {*************************  soper  ********************}
      SUM := 0;
      I := 0;
      AI := 1;
      V := X /XL;
      R := B + XL * (A + B);
      SMAX := MAXINT - 1;
      If R < SMAX Then
         S := TRUNC(R)
      Else S := SMAX;
      If S > 0 Then
         Repeat
         {****************************************************}
         {**    Raise a Decreas b                           **}
         {****************************************************}
            SUM := SUM + AI;
            I := I + 1;
            AI := AI * V * (B - I ) /(A + I);
         Until ( Abs(AI /SUM) < EPS) or ( I >= S);
      If ( I = S ) Then
         Begin
         {**************************************************}
         {**    Raise a                                   **}
         {**************************************************}
             K := 0;
             CK := XL * AI;
             ABL := A + B - 1;
             AS := A + S;
             Repeat
                SUM := SUM + CK;
                K := K + 1;
                CK := X * CK * (ABL + K) /(AS + K);
            Until (ABS(CK /SUM) < EPS);
      End; { raise a}
      A12 := (A + 1) * (A + 2);
      B2 := B * (B + 1);
      AB4 := (A + B) * (A + B + 1) * (A + B + 2) * (A + B + 3);
      SOPER := SUM * F(X, XL, A + 2, B + 2) * A12 *
          B2 /AB4 /SQR(X * XL) /XL;
   End; {soper}

Begin
{****************************  ICBETA  *******************}
If X <= 0 Then
   ICBETA := 0
Else
   If X >= 1.0 Then
      ICBETA := 1.0
   Else
      If (A > 1) And (B > 1) Then
         If X <= (A /(A - B)) Then
            ICBETA := SOPER(X, XL, A, B)
         Else ICBETA := 1.0 - SOPER(XL, X, B, A)
      Else
         Begin
         {********************************************}
         {**    Increment a,b                       **}
         {********************************************}
         C:= A * (A + 1) /(A + B) * (A + 2) /( A + B + 1) *
            B /(A + B + 2) * (B + 1) /(A + B + 3);
         FXAB := F(X, XL, A + 2, B + 2) * C /Sqr(X * XL);
         ICBETA := FXAB * (1 - X * (A + B) /B) /A +
            ICBETA(X, XL, A + 1, B + 1);
         End;
      {**********************************************}
      {**     Increment a,b                        **}
      {**********************************************}
End;  { ICBETA}

Begin
{***************************  Probability  ***************}
      N := TRIALS;
      P  := 0.2;
      PL := 0.8;
      X := HITS;
      PROB_B := ICBETA(PL, P, N - X, X + 1);
End;          {probabilty}

Procedure CHISQUARE(Var PROB_C : Real;
                        CHI_SQ : Real;
                        DEGREE_FREEDOM : Real);
Var
X,N : real;
A : Real;
TEMP : String[18];
{*********************************************************}
{**                                                     **}
{**     Procedure CHISQUARE                             **}
{**                                                     **}
{*********************************************************}
Function BINET(A:real) :Real; { Benet's Function }
   Var
   S : Real;
   Begin {BINET}
   {*************************  binet  ********************}
      If A < 1.5 Then
         BINET := 8.1061467E-2
        {*************************************************}
        {**      Case BINET(1) = 1.5*ln(2pi)            **}
        {*************************************************}
     Else
        {*************************************************}
        {**    Continued Fraction                       **}
        {*************************************************}
        Begin
           S := A + 1.5174736 /(A + 2.269489);
           S := A + (195 /371) /(A + 22999 /22737 /S);
           BINET := 1 /12 /(A + (1 /30) /(A + (52 /210) /S));
        End;
End; {BINET}

Function G(U : Real) : Real;
{*********************************************************}
{**   ln(1+u)-u with accuracy                           **}
{*********************************************************}
Var         R  : Real;
Begin {G}
   {  book corrected)
   If (U < -0.5) Or (U > 1.0) Then
      G := Ln(1 + U) - U
   Else
   {*************************************************}
   {**   continued fraction                        **}
   {*************************************************}
   Begin
      R := 6 + 16 * U /(7 + 9 * U /(8 + 25 * U /(9 + 16 * U /10)));
      G := -U * U /(2 + 4 * U /(3 + U /(4 + 9 * U /(5 + 4 * U /R))));
   End;
End;  { G}

{************************************************}
{**                                            **}
{**    Function LNGAM2                         **}
{**                                            **}
{************************************************}
Function LNGAM2(A:Real) : Real ; { LN(GAMMA(2 + A)) }
Var
K : Integer;
S : Real ;
C : Array[1..9] Of Real; { C[1]=1-GAMMA, C[K]=(-1**K)*(ZETA(k)-1)/k}
{******************  LNGAM2  *******************}
Begin {LNGAM2}
   C[1] :=  0.4227843; C[2] :=  0.3224670; C[3] := -0.0673523;
   C[4] := -0.0205808; C[5] := -0.0073856; C[6] :=  0.0028905;
   C[7] := -0.0011928; C[8] :=  0.0005097; C[9] := -0.0002232;
   S := C[9] * A;
   For K := 8 DownTo 1 Do S:= C[K] + A * S;
      LNGAM2 := A * S;
End; {LNGAM2}

{************************************************}
{**                                            **}
{**    Function ICGAMMA                        **}
{**                                            **}
{************************************************}
Function ICGAMMA(A,X : Real) : Real; {Incomplete GAMMA C.D.F.}
Var
K: Integer;
D,S,C,UK,AK : Real;
{************************************************}
{**                                            **}
{**    Function LNXF                           **}
{**                                            **}
{************************************************}
   Function LNXF(X,A : Real) : Real; {Ln(X*F(X,A))}
   Var
   U,GU : Real;
   {********************  LNXF   ****************}
   Begin {LNXF}
      U := (X - A) /A;
      If U < -0.9 Then GU := Ln(X /A) - U
      Else GU := G(U);
      If A < 0.5 Then LNXF := A * GU + G(A) + (A + 1) * LN(A) - LNGAM2(A)
      Else
         If A<1.5 Then
            LNXF := A * GU + G(A - 1) + (A * LN(A) - 1) - LNGAM2(A - 1)
         Else {case a >= 1.5}
            LNXF := A * GU - BINET(A) - 0.5 * Ln(6.2831853 /A);
   End; {LNXF}

{********************  ICGAMMA  **************}
Begin {ICGAMMA}
   If X <= 0 Then ICGAMMA := 0
   Else
      Begin {nonzero case}
         D := Exp(LNXF(X,A));
         If (X <= A) Or (X <= 1) Then
            Begin {series}
               AK := A;
               UK := 1;
               S := 1;
               While (Abs(UK /S) > 1E-7) Do
                  Begin
                     AK := AK +1;
                     UK :=UK * X /AK;
                     S := S + UK;
                  End;
               ICGAMMA := S * D /A;
            End {series}
         Else
            Begin {continued fraction}
               C:= 0;
               For K := 30 DownTo 1 Do
                  Begin
                     C := K /(X + C);
                     C := (K - A) /(1 + C);
                  End;
               C:= 1 /(X + C);
               ICGAMMA := 1 - C * D;
            End;  {continued fraction}
      End;  {non zero case}
End; {ICGAMMA}

Begin { CHISQUARE}
{**************************  CHISQUARE  *****************}
N :=  DEGREE_FREEDOM;
X :=  CHI_SQ;
If N < 0.5 Then Halt(0)
Else
   Begin
      If X < 0 Then PROB_C := 0
      Else PROB_C := ICGAMMA(N /2.0,X /2.0);
   End;
End; { CHISQUARE }

End.    { Turbo Pascal Unit  --  ESP_PROB}